{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import SimpleITK as sitk\n",
    "print(sitk.Version())\n",
    "from myshow import myshow\n",
    "# Download data to work on\n",
    "%run update_path_to_download_script\n",
    "from downloaddata import fetch_data as fdata\n",
    "\n",
    "OUTPUT_DIR = \"Output\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This section of the Visible Human Male is about 1.5GB. To expedite processing and registration we crop the region of interest, and reduce the resolution. Take note that the physical space is maintained through these operations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fixed_rgb = sitk.ReadImage(fdata(\"vm_head_rgb.mha\"))\n",
    "fixed_rgb = fixed_rgb[735:1330,204:975,:]\n",
    "fixed_rgb = sitk.BinShrink(fixed_rgb,[3,3,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "moving = sitk.ReadImage(fdata(\"vm_head_mri.mha\"))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myshow(moving)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Segment blue ice\n",
    "seeds = [[10,10,10]]\n",
    "fixed_mask = sitk.VectorConfidenceConnected(fixed_rgb, seedList=seeds, initialNeighborhoodRadius=5, numberOfIterations=4, multiplier=8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Invert the segment and choose largest component\n",
    "fixed_mask = sitk.RelabelComponent(sitk.ConnectedComponent(fixed_mask==0))==1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myshow(sitk.Mask(fixed_rgb, fixed_mask));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# pick red channel\n",
    "fixed = sitk.VectorIndexSelectionCast(fixed_rgb,0)\n",
    "\n",
    "fixed = sitk.Cast(fixed,sitk.sitkFloat32)\n",
    "moving = sitk.Cast(moving,sitk.sitkFloat32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "initialTransform = sitk.Euler3DTransform()\n",
    "initialTransform = sitk.CenteredTransformInitializer(sitk.Cast(fixed_mask,moving.GetPixelID()), moving, initialTransform, sitk.CenteredTransformInitializerFilter.MOMENTS)\n",
    "print(initialTransform)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def command_iteration(method) :\n",
    "    print(\"{0} = {1} : {2}\".format(method.GetOptimizerIteration(),\n",
    "                                   method.GetMetricValue(),\n",
    "                                   method.GetOptimizerPosition()),\n",
    "              end='\\n');\n",
    "    sys.stdout.flush();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tx = initialTransform\n",
    "R = sitk.ImageRegistrationMethod()\n",
    "R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)\n",
    "R.SetOptimizerAsGradientDescentLineSearch(learningRate=1,numberOfIterations=100)\n",
    "R.SetOptimizerScalesFromIndexShift()\n",
    "R.SetShrinkFactorsPerLevel([4,2,1])\n",
    "R.SetSmoothingSigmasPerLevel([8,4,2])\n",
    "R.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()\n",
    "R.SetMetricSamplingStrategy(R.RANDOM)\n",
    "R.SetMetricSamplingPercentage(0.1)\n",
    "R.SetInitialTransform(tx)\n",
    "R.SetInterpolator(sitk.sitkLinear)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import sys\n",
    "R.RemoveAllCommands()\n",
    "R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )\n",
    "outTx = R.Execute(sitk.Cast(fixed,sitk.sitkFloat32), sitk.Cast(moving,sitk.sitkFloat32))\n",
    "\n",
    "print(\"-------\")\n",
    "print(tx)\n",
    "print(\"Optimizer stop condition: {0}\".format(R.GetOptimizerStopConditionDescription()))\n",
    "print(\" Iteration: {0}\".format(R.GetOptimizerIteration()))\n",
    "print(\" Metric value: {0}\".format(R.GetMetricValue()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tx.AddTransform(sitk.Transform(3,sitk.sitkAffine))\n",
    "\n",
    "R.SetOptimizerAsGradientDescentLineSearch(learningRate=1,numberOfIterations=100)\n",
    "R.SetOptimizerScalesFromIndexShift()\n",
    "R.SetShrinkFactorsPerLevel([2,1])\n",
    "R.SetSmoothingSigmasPerLevel([4,1])\n",
    "R.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()\n",
    "R.SetInitialTransform(tx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "outTx = R.Execute(sitk.Cast(fixed,sitk.sitkFloat32), sitk.Cast(moving,sitk.sitkFloat32))\n",
    "R.GetOptimizerStopConditionDescription()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "resample = sitk.ResampleImageFilter()\n",
    "resample.SetReferenceImage(fixed_rgb)\n",
    "resample.SetInterpolator(sitk.sitkBSpline)\n",
    "resample.SetTransform(outTx)\n",
    "resample.AddCommand(sitk.sitkProgressEvent, lambda: print(\"\\rProgress: {0:03.1f}%...\".format(100*resample.GetProgress()),end=''))\n",
    "resample.AddCommand(sitk.sitkProgressEvent, lambda: sys.stdout.flush())\n",
    "resample.AddCommand(sitk.sitkEndEvent, lambda: print(\"Done\"))\n",
    "out = resample.Execute(moving)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "out_rgb = sitk.Cast( sitk.Compose( [sitk.RescaleIntensity(out)]*3), sitk.sitkVectorUInt8)\n",
    "vis_xy = sitk.CheckerBoard(fixed_rgb, out_rgb, checkerPattern=[8,8,1])\n",
    "vis_xz = sitk.CheckerBoard(fixed_rgb, out_rgb, checkerPattern=[8,1,8])\n",
    "vis_xz = sitk.PermuteAxes(vis_xz, [0,2,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "myshow(vis_xz,dpi=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "sitk.WriteImage(out, os.path.join(OUTPUT_DIR, \"example_registration.mha\"))\n",
    "sitk.WriteImage(vis_xy, os.path.join(OUTPUT_DIR, \"example_registration_xy.mha\"))\n",
    "sitk.WriteImage(vis_xz, os.path.join(OUTPUT_DIR, \"example_registration_xz.mha\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
